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ABSTRACT 


A new finite element modeling is presented to 
investigate the static and dynamic behavior of laminated 
composite beams with partial delamination. In this study, a 
newly developed rectangular beam element is used. The element 
has lateral and axial displacements as degrees of freedom but 
no rotation. For simplicity, linear shape functions are used 
for the beam element. As a result, the element has six degrees 
of freedom, four of which are the axial displacements at the 
corner points and two are the lateral displacements at the 
ends. In addition, contact-impact conditions are applied to 
the finite element modeling to avoid overlapping of the upper 
and lower portions of a delaminated section. The numerical 
study shows that depending on existence of an embedded 
delamination crack and its size, the response is different for 


a beam with a crack and subjected to a short impulse load. 
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I. INTRODUCTION 


Laminated composite structures provide improved mechanical 
properties such as high strength, stiffness, toughness and 
high temperature performance [Ref . Le) However, these 
structures have more complicated failure mechanisms and modes 
compared to those of conventional metallic structures. One of 
the common failure mode is inter-delamination. 

The objectives of this study are to develop a finite 
element model for an analysis of laminated composite beam 
structures containing an inter-delamination crack and to find 
a methodology to detect an embedded interlaminar delamination. 
For delaminated beams, a short impulse load is applied to find 
out whether the dynamic response varies significantly compared 
to the dynamic response of structures without a crack. 

In this study, a newly developed finite element model 
[Ref. 2] is used for the simulation of beams containing an 
inter-delamination crack. A rectangular element that has six 
degrees of freedom is used. The element allows axial 
displacements at its four corner points and lateral 
displacements at its two ends. Using this element, static and 
dynamic analyses of laminated composite beams containing a 
local delamination crack are conducted. First of all, static 
deflections as well as natural frequencies and mode shapes of 


beams without delamination are calculated to check the 


accuracy of the element. Laminated composite beams with a 
partial delamination are studied using the contact-impact 


conditions applied at the delaminated portion. 


II. FINITE ELEMENT FORMULATION 


A. FINITE ELEMENT MODELING 

A newly developed rectangular finite element, Figure 1, is 
used to investigate the static and dynamic behavior of 
laminated composite beams. The element has six degrees of 
freedom and allows axial and lateral displacements but no 
rotation. There are axial displacements at the four corner 


points and lateral displacements at the two ends of the 








element. 
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Figure 1 :Four Noded Rectangular Element with Six Degrees of Freedom. 


In Figure 1 and the subsequent formulation, u represents 


the axial displacements at the corner points and v represents 
the lateral displacements at the ends. The subscripts "1" and 
"2" refer to the left and right ends while the superscripts 
"t" and "b" indicate the top and bottom sides of the element, 
respectively. 


The displacement model of the element is 


v(x) 


a-{250)}- cide) (1) 


where [N] is the matrix of shape functions and {d°*} is a 
vector of nodal displacements. The axial displacement is 
assumed to vary linearly along both axial and lateral 


directions. It can be written as 


2 

u(x, Y)= 2 Ni (x) (LH, (y) u?+H,(y) u,'] (2) 
= WN, (2) Hy (y) uP +N, (x) H, (y) uy +N, (x) H, (y) 0? +N, (x) H (y) u," 

The lateral displacement, which is assumed to be constant 


through the thickness of the element, varies linearly along 


the axial direction and can be written as 


2 
S2N; : 
v(x) i (x) V; (3) 
=N, (x) v,+N, (x) v, 


Here, N, and H,; are the linear interpolation or shape functions 
in the axial and lateral directions. The beam element may use 
a higher order shape function for the axial direction, i.e. 
N,, and the linear shape function for the lateral direction, 


i.e. H,. However, in this study, the linear shape function is 


used for both N, and Hi for simplicity. That is, 


N, (x) eres 
a 


N, (x) == 
% (4) 
B,(y) =1-4 
a a 
H,(y) aay 


where a and b are the length and height of the beam element, 
respectively. 

For simplicity, notations N,, N,, H, and H, will be used 
instead of N,(x),N,(x),H,(y) and H,(y) in the following 
derivation. 

For two dimensional elasticity problems the strain vector 


can be expressed as [Ref. 3] 


e,| la/ax 0 
B=}€, (=| 0 ad/ay|r= [Bas (5) 
xy} |8/ay a/dx 


where {d°} is the vector of nodal displacements and [B] is the 
matrix that relates strains to the nodal displacement vector. 


Normal strains can be written as 


gu OM 4 Ou Om, OM 

Pe Bae Oa BY * Gag MaMa * Gag Tha * Gi Halt 

es ov_ (6) 
i oy 


and the shear strain is 








The element stiffness matrix can be obtained by minimizing 
the total strain energy which contains both bending and the 
transverse shear energy. This minimization yields the 


following element stiffness matrix 
[k*]=[k?]+[k?] (8) 


where the subscripts "B" and "S" indicates bending and the 
transverse shear respectively. 
The bending and transverse shear stiffness matrices are 


given as 


ab 

[k= f fi (B) E(B, dydx (9) 
00 
ab 

lke }= f f {(B}7G(B) dydx (10) 
00 


where E and G are the elastic and shear modulii of the beam 
and the vectors (B,} and (B,} are derived below. 
The strain-displacement relationship for the normal strain 


in the axial direction is 


OH, 4, OH, aM 


(11) 


(12) 


(13) 


The bending stiffness matrix can be obtained by carrying 


(9) which will result in 


out the integration in Eq. 








For Eq. (10) the reduced integration technique is used 
along the x-axis to prevent shear locking which occurs when 


the ratio of beam length to beam thickness is large. This 


integration yields the transverse shear stiffness matrix of 














the form 

[Ga Ga G Ga Ga Gl 

4b 4b 206 4BOC BD 

_Ga Ga _G _Ga Ga G 

4b 4b 2 4b Ab 2 

G G@ ¢b & G Gb 

aoe 2 2 a 2 ar 

sl") Ga Ga G Ga Ga ue 

4b 4 2 4b 4b 2 

_Ga Ga _G Ga Ga G 

4b 4b 2 4b 4b 2 

_& G _Gbh _G@ G Gb 
| 2 2 a 2 8 a] 


The element stiffness matrix, which is obtained by adding 





the bending and transverse stiffness matrices, can be 
expressed in the following form 
a,+2a, -aj+a, a, a,-2a, -a,-a, -a 
-a,+@, a+2a, -a, -a,-a, a,-a, a 
a -a a -a -a 
4 4 2 4 4 2 
[k*]= (16) 
a,~2a, -a,-a, a a)+2a, -a,ta, —a, 
“@,-a, a,-2a, -a, -a,ta, ata, a, 
“ay a, “a, “ay a, a 
where 
a = 28> y-FR 4.6 
ap 92°] 3 6a + 2 (17) 


The mass matrix can be derived consistently in a similar 
way. In this study a diagonal (lumped) mass matrix was used 


since it takes less computing time. The mass matrix is 


where p is the mass density, a and b are the length and height 


of the element. The element is assumed to have a unit depth. 


B. LAMINATED COMPOSITE BEAMS 


In the most generalized form Hooke's law can be written as 


O55=Cina® ia tad pke t= 18,5 (19) 


using tensor notations [Ref. 4]. Here C,;,, represents the 
elastic stiffness. Using the symmetry of stress and strain 
tensors, i.e. 6,,-0,; and EE) and using a shorthand notation 


we can write the normal stresses and strains as 


O1=011 Ey=Far 





O2=022 E2=F22 
03=033 E,=F33 


while the shear stresses and strains are written as 


O4=023-032 E4g=F 232832 
O5=031-913 Es=E 315813 
Og=9) 2-921 E6=Ep22F 21 


For general orthotropic materials Hooke's law can be written, 


in the inverted form, as 


[s,, Sin Siz 9 +O 0 
€ iO 
- Sy Sy, 0 0 0 If 
2 2 
esl. S, 0 0 0 /|Io, (0) 
e 4) Si, 9 0 (1a 
& oO 
5 ic 0 5 
&, 55 O~ 
Se6| 








where [S,,] is the compliance matrix, which is the inverse of 
the matrix [C,,]. For the state of plain stress, Eq. (20) is 


reduced to 

‘ Si, Siz 0 2 

Eof=|S12 So. 0 |9o, (21) 
J |0 0 s,,|C 

Stress-strain relation can be obtained by inverting the 


relation given in Eq. (21) 


co, Q, 2, 0 S. 
S2f=|Q2 OQ. O |)€2 (22) 


J i/o 0 Q,|Fs 
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where the reduced stiffness terms are 


Qa 
ee 2 
Wr aa 
1Y2 (23) 
_ VE, VE, 
RB d-v,v, 1-v,), 
2667 Se 


In this study, the width of the beam is neglected, i.e. it is 
assumed to have a unit width. Therefore, only the elastic 
modulus along the beam axis is considered. Two different 
laminae with four and eight layers are used for calculations, 
l.e. [0/90], and [0/90/0/90], where the subscript '"s" 
represents symmetry with respect to the middle axis of the 
beam. 

There are several methods used in modeling laminated 
composite beams. In this study, two different techniques are 
used. The first technique discretizes respective layers in the 
finite element analysis. As a result, the total number of 
elements is proportional to the number of layers in a 
laminated beam. Of course, static condensation can be 
performed to reduce the number of total degrees of freedom by 
eliminating internal layers degrees of freedom. This modeling 
technique is computationally expensive but it can describe a 
general shape of deformation through the beam thickness. The 
second technique uses one beam element through the beam 


thickness regardless of the number of layers. This technique 


1] 


is computationally efficient but it assumes aé_ linear 
deformation through the beam thickness. The development of 


this technique is described in the following paragraphs. 























F igure 2 : Laminated Beam Element 


Let u, and v, represent the axial and lateral displacements 
of the kth layer and let h be the beam thickness while h, and 
h,,, represent heights of the top and bottom sides of the layer 
measured from the bottom of the beam (see Figure 2). The 
relationship between the layer displacements and the global 


beam displacements can be written as 


b _ b t 
Wy Guy Foy 
at <6 wheat 
bea Cay FOg Ry 








Yea My 
gh =¢, ou" 04) 
yt b 
Ux. 2 = ¢,U, Heyy 
Vi,2~ V2 
where the constants, c's, are 
h-h 
c\= ss 
h 
h 
or 
h-h (25) 
= kel 
or 
h 
c= ke 
A 
This relationship can be expressed in the matrix form as 
Go 0 0 0 
ge ~ 38 uy” 
gs c, Cc, 0 0 0 uf 
V,.4{ |0 O 1 0 O}}v, 
B f= 26 
ms) (OO Oe ce 0 le 28) 
Uso! |0 0 0 c, c, 0] 
V, V. 
e2} 19 0 0 0 O 1)? 
or 
{d= [T]{d} (27) 


where {d*} and {d} are displacement vectors of the kth layer 
and the beam, and [T] is the transformation matrix. Now, 


stiffness and mass matrices of a laminated beam element can be 


expressed as 


13 


key =3° [TIT(k*] (T] 
k=1 


(28) 
[m°] n% [T]7{m*] [7] 
k=1 


Here, [k*] and [m‘] are the stiffness and mass matrices of the 
kth layer, respectively. Using this technique, a single beam 


element can include all the layers of a laminated beam. 


C. DYNAMIC ANALYSIS 

In dynamic analyses, displacements, velocities, 
accelerations, stresses and strains need to be calculated at 
each time interval since they all vary with time. This 
requires a time integration scheme. In this study, a form of 
central difference method is used in the dynamic analysis of 
undamped systems. 


The general equation of motion at time t can be written as 


[M] {u} + [K]{u},={F), (29) 


The central difference method can be applied in three steps. 
A) Compute [M] and [K] 


B) Initial calculations 


1. {u),=[™M] + (F},- [K]{u},) 


ee {a} ae={u) ae oui}, (A) (30) 


3. {uly ={u)ad ar (At) 
2 


C) For each time step 


1. {u),=[M] > ({F) ,- [K]{u),) 
ai fa), actu, oe dud} (Ae) (31) 


3. ful... =tub+{a), oc (At) 
2 


This method is conditionally stable when the time step 
Size At is less than or equal to the critical step size [Ref. 


5,6] 


a= (32) 


where (@,,)° is the largest eigenvalue of the matrix equation 
(19) with no forcing function. When At>At,., it results in an 


unstable numerical solution. 


D. CONTACT-IMPACT PROBLEMS 

In the case of partially delaminated beams, as shown in 
Figure 3, the upper and lower portions of the delaminated 
section may overlap each other if there is no constraint 
applied on the portion during the finite element calculations. 
This fact applies to both static and dynamic problems. To 


avoid this overlap the contact-impact conditions are applied. 
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Figure 3 : Simply Supported Beam with a Delaminated Section. 


The contact-impact algorithm developed in [Ref. 7] is 
applied to the present beam element as explained in the 
following paragraphs. 

Let the letters i and j represent the facing nodal points 
in the upper and lower portions of a delaminated section, and 
let v, and v, represent the vertical displacements of these 
points. The procedure to apply the contact-impact condition is 
explained below. 

A) Initially, the central difference method is used to 
calculate displacements, velocities and accelerations until 
any two points contact each other (v,-v,>0). 


B) When the two nodes penetrate each other, i.e. v,-v,<0, the 
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release-to-contact condition is applied. And the corrected 
values of the velocity, acceleration and elastic force are 


calculated using the following formula; 


L6 A. 2542 
ee u_,tM uly 


uw =—___— (33) 





C) When they are released from each other, i.e. V;-V;>0, the 
contact-to-release condition is applied and the corrected 
values of velocity and acceleration are computed using the 


following formula; 


(=L)"ACt., 
i ry 
(-1)@1_ (34) 
ways _ 
we mM 


where the value of a is 1 for the upper portion and 2 for the 
lower portion of the delaminated section. Here the subscript 
(+) is used for the updated values while the subscript (-1) 
indicates the values taken at the end of the previous time 
step and (-) indicates the values taken at the end of the last 
iteration of the present time step. 

For static problems, the contact condition is applied 
simply by setting the vertical displacement of the upper 
portion equal to the vertical displacement of the lower 


portion where they overlap each other. 
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III. RESULTS AND DISCUSSION 


A. VERIFICATION EXAMPLES 

The present finite element model and the computer code are 
verified by analyzing three problems and comparing the results 
with exact solutions. The verification examples are static 
deflections of both isotropic and laminated composite beams 
and the free vibration of beams. Two different boundary 
conditions are applied. One of them is a simply supported beam 


and the other is a cantilever beam. When loading is required, 








a) Pp 


/\ ( ) 
haan Z WV 
b) P 














Figure 4: a) Simply Supported Beam with Center Load b) Cantilever Beam with End 
Load. 
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the force is applied at the center of the simply supported 
beam and at the end of the cantilever beam as seen in Figure 
4. 

Material properties of the beams are as follows; elastic 
modulus E=1,000,000 Pa., shear modulus G=450,000 Pa., density 
o=10 g/cm’, length of the beam L=60 m., height h=0.5 m. and 
applied force P=1 N. 

First a free vibration analysis is performed, for both 
simply supported and cantilever beams. The analytical 


expression for the natural frequencies is as follows 


>. (k,1)4ET 
O.=—$—$____—- 


ya (35) 
where y is the mass per unit length of the beam and k,l are 
defined values for some beam types. 

The results of first three natural frequencies are 
presented in Tables 1 and 2. As seen from the results, the 
finite element solution agrees with the exact solution very 
well. The mode shapes corresponding to first four natural 
frequencies are also plotted in Figures 5 and 6. These are 
expected mode shapes of the simply supported and the 


cantilever beams [Ref. 8]. 
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TABLE 1 : NATURAL FREQUENCIES OF SIMPLY SUPPORTED BEAM 




















TABLE 2 : NATURAL FREQUENCIES OF CANTILEVER BEAM 











Then static deflections of beams are computed. The 


analytical solutions for isotropic beams are [Ref. 9] 








mS -PxX 2 2 L 
-Px? 
ve (Bix) (37) 


where Eq. (36) is for a simply supported beam and Eq. (37) is 
for a cantilever beam. Here "I" refers to the moment of 
inertia and L is the length of the beam. The results of both 


exact solution and finite element solution are plotted in 


Figure 7. 
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Figure 5 : Mode Shapes of the Simply Supported Beam 
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Figure 6 : Mode Shapes of the Cantilever Beam 
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Figure 7 : Static Deflections of Isotropic Beams; a) Simply Supported Beam , 
b)Cantilever Beam 
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Finally a static analysis of a laminated composite beam 
with simple supports is performed. The analytical expression 
for the deflection of a simply supported beam, which is 
subject to three point bending, as shown in Figure 8, is as 


follows [Ref. 10] 








pe { 2) ° 38 
48EPT L an 
and the equivalent bending stiffness is defined as 
Bee geek 
EPI=D ESI (39) 


k-1 


here E,° is the effective bending modulus of the beam, E,* is 
the modulus of the kth layer relative to the beam axis, I“ is 
the moment of inertia of the kth layer relative to the 
midplane, and N is the number of layers in the laminate. The 
vertical displacements obtained from the exact solution and 
the finite element solution are plotted together, for four and 
eight layered beams and for two different ratios of elastic 


modulii, in Figures 9 and 10. 
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Figure 8 : Laminated Beam with Three Point Bending (Simply Supported). 

In the static problems, as seen from the figures, the 
exact solutions and the finite element solutions agree very 
well. From these verification examples, it can be shown that 
the beam element provides accurate results compared to exact 


solutions. 
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Figure 9 - Static Deflections of the Laminated Composite Beam with 4 Layers 
a)for E,/E,)=20, b) for E,/E,,=100 
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Figure 10: Static Deflections of the Laminated Composite Beam with 8 Layers 
a) for Ep/E,9=20, b) for E,/E,>=100 
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B. NUMERICAL RESULTS 

Dynamic analyses are performed for a simply supported 
composite beam with four layers, i.e. [0/90],, and with 
different crack sizes using contact-impact conditions. 
Calculations without contact-impact conditions are also made 
and compared to the case computed with contact-impact 
conditions. An impulsive load of a very short duration is 
applied at the center of the upper portion of a delaminated 
section and the displacement response of the beam is recorded 
for a period of time. 

First displacements at the mid-points of upper and lower 
portions of a large crack are calculated and plotted in Figure 
11 without using contact-impact conditions, where crack length 
is one third of the beam length. As seen from the figure, 
these two points overlap each other which is not admissible. 
Then for the same crack size, contact-impact conditions are 
applied and displacements of the beam is calculated and 
plotted. Figure 12 shows the displacements of the beam just 
before contact while the displacements at the time of contact 
is shown in Figure 13. The displacements just after contact is 
shown in Figure 14. Also the responses of the mid-points and 
quarter-points are plotted in Figures 15 and 16, respectively, 
where quarter-points represent the points which are located in 


the middle between the crack tip and the crack midpoint for 
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both upper and lower portions. Later, the crack size is 
reduced to one thirtieth of the beam length and the 
displacements of mid-points and quarter-points are plotted in 
Figures 17 and 18, respectively. As seen in these results, the 
upper and lower portions never penetrate each other. This 
indicates that the application of contact-impact conditions 


works well in simulating delaminated beams. Finally, the 


response of a beam without a crack is computed and plotted in 


Figure 19 for comparison purpose. Responses of all these cases 
are different each other. They may give an idea how an 


embedded crack affects the behavior of beams. 
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Figure 11 : Displacements of Mid-Points without Contact-Impact Conditions 
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Figure 12 | Displacement of Beam (before Contact) 
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Figure 13 : Displacement of Beam (at Contact) 
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Figure 14 : Displacement of Beam (after Contact) 


33 








-4 
x 10 





<== upper portion 


displacement, m. 


<== lower portion 








time, sec. 











Figure 15 : Displacements of Mid-Points with Contact-Impact Conditions 
(Large Crack Size) 
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Figure 16 : Displacements of Quarter-Points with Contact-Impact Conditions 
(Large Crack Size) 
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Figure 17 : Displacements of Mid-Points with Contact-Impact Conditions 
(Small Crack Size) 
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Figure 18 : Displacements of Quarter-Points with Contact-Impact Conditions 
(Small Crack Size) 
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Iv. CONCLUSIONS AND RECOMMENDATIONS 


The beam element used in this study is useful to model a 
crack with contact-impact conditions. Application of contact- 
impact conditions constrains the vertical displacements at the 
crack in a way that the points on lower and upper portions of 
the crack hit each other without penetration. 

Computational time depends on the number of degrees of 
freedom. The technique, which models the laminated composite 
beam using a single beam element, reduces the degrees of 
freedom to save time during calculations. 

When an impulsive load of a short duration is applied, the 
beam responds in different ways depending on whether there is 
a crack inside it and what the crack size is. Examining these 
results gives an idea about existence of an embedded crack as 
well its size and location. 

For further studies, delaminated sections may be refined 
to see the minimum crack size that can be detected using this 
procedure. Also, the Fast Fourier Transform can be applied to 
measure the frequency response of a delaminated beam. 
Comparing the natural frequencies may give an idea about the 
crack size. In this study, a crack is assumed to be along the 
middle axis of the beam and the load is applied at the center 
of the beam. Locations of cracks and loads may be changed and 


the responses can be calculated for different cases. 
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Some experiments are being performed to see the response 
of beams under impulsive loading. The results obtained from 
the procedure described in this study may be compared with 


those experimental results for verification. 
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